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L8: Entry 1 of 6 File: USPT Oct' 23, 2001 

DOCUMENT -IDENTIFIER: US 6306092 Bl 

TITLE: Method and apparatus for calibrating rotational offsets in ultrasound 
transducer scans 



ABPL: 

The axis of rotational transducer array grans, because of imperfect transducer 
array assembly, may have two orthogonal offsets relative to the geometric 
center of the transducer array. Without knowledge of these offsets, it is not 
possible to convert rotational transducer scan data into a rectilinear 
(Euclidean) coordinate system, as is necessary for three-dimensional 
processing. Using spatial coherency between appropriate scan lines in 
different rotational transducer scans, the horizontal and vertical rotational 
offsets are calculated. These offsets are then utilized in converting the data 
to a rectilinear coordinate system for three-dimensional processing. 

BSPR: 

Three-dimensional visualization has gained popularity in medical applications 
since the introduction of computer tomography (CT) in the field several 
decades ago. For example, three-dimensional visualization is also used in 
magnetic resonance. (MR) imaging. Using three-dimensional data sets in 
ultrasound imaging is not as popular due to two major obstacles: first, data 
in most cases are acquired by free-hand B-mode scans that do not provide 
sufficiently accurate information to enable precise positioning of individual 
two-dimensional scans (slices) into a common three-dimensional coordinate 
space; second, the ultrasound data are inherently more noisy than CT and MR 
data sets, and therefore traditional' surface visualization techniques do not 
produce good results. The last decade has brought many advances in technology, 
in both hardware and software, that allow for real-time three-dimensional data 
set visualization using so-called volume rendering that goes directly from a 
three-dimensional data set into a two-dimensional image display, bypassing the 
creation of surfaces. One volume rendering technique is known as maximum 
intensity projection (MIP) . The MIP technique involves projection of 
three-dimensional data intensity values onto a two-dimensional image plane by 
assigning to each image pixel the maximum intensity value in the 
three-dimensional data volume that belongs to the line of sight that goes from 
the eye point through this pixel and into the volume. This method, in 
combination with animation, can produce true three-dimensional impressions on 
the monitor. A more computationally demanding technique is known as 
compositing. This technique involves modeling of a physical phenomenon of 
light propagation in semi -trans lucent /semi -opaque media that is recreated from 
a three-dimensional data set with the addition of specially designed transfer 
functions . 

BSPR: 

Some medical applications involve acquiring three-dimensional volume data by a 
transducer that rotates around an axis orthogonal to the transducer array. The 
volume "swept" by these two- dimensional B-mode scans represents a cylinder. 
Since the two-dimensional scans do not lie parallel to each other, it is 
difficult to visualize three-dimensional object structures from individual 
scans alone and a volume visualization technique would be desirable. 
Contemporary software and hardware are efficient in volume rendering 
techniques, but require that the data be represented as a rectilinear 
three-dimensional data array. Therefore, conversion from a cylindrical 
coordinate system to a rectilinear coordinate system is required. Although 
this conversion is not difficult to compute, an important practical 
complication to the conversion process is that there is always some offset of 
the axis of rotation relative to the sensor array middle point. Need exists 
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c 



for a method of calcTFlating this offset based on one^nree- dimensional volume 
scan. The offset information computed can be used in an algorithm to convert 
from a cylindrical coordinate system to a rectilinear one and also can be used 
in the transducer manufacturing process to position a transducer array exactly 
at the rotational axis. 



DRPR: 

FIG. 9 is a schematic illustration of two overlapping coordinate grids in a 
plane orthogonal to the rotational transducer scans ; one induced by rotational 
scans and the other induced by parallel scans . 

DRPR: 

FIGS. 17 and 18 are schematic representations of two consecutive rotational 
transducer scans . 



DRPR: 

FIG. 19 is a graph of the correlation function between corresponding vertical 
scan lines in two consecutive rotational transducer srans . 



DRPR: 

FIGS. 22 and 23 are schematic representations of two. rotational transducer 
scans corresponding to slices PQ and P'Q', respectively, in FIG. 21. 

DRPR: 

FIG. 24 is a graph of the correlation function between vertical sran lines in 
two conjugate slices (PQ and P'Q'). 

DEPR: 

Rotational transducer scans are generated by a one -dimensional transducer 
array mounted on a rotating plate 3 (FIG. 1) having an axis of rotation 
perpendicular to the transducer array. The invention addresses the situation 
wherein the center of the transducer array is offset from the axis of rotation 
of the rotating plate. Each aran represents a two -dimensional data array 36, 
shown in FIG. 4, which is M elements wide and N elements high. The totality of 
scans for all rotational angles of the transducer sweeps a cylindrical volume 
38, shown in FIG. 3. Various post-processing techniques usually deal with 
rectilinear structured grid sets. FIGS. 5 and 6, respectively, represent a 
three-dimensional rectilinear volume data set 40 and one of its 
two-dimensional horizontal slices 42. The slices are M elements wide and M 
elements deep. The volume is N elements high. The total amount of elements in 
the volume data set 4 0 is M. times .M. times .N . When converting from rotational 
slices to rectilinear coordinates, sampling over volume M . times .M. times . N need 
not be performed with the same sampling steps along rectilinear axes. However, 
without loss of generality, the same number of elements along horizontal 
rectilinear axes will be used as the number of elements in the rotational 
transducer. 



DEPR: 

In accordance with a preferred embodiment of the invention, the procedure for 
computing offset_x is as follows. In FIG. 15, the lines PQ and P'Q' represent 
top views of two consecutive rotational transducer scans FIG. 16 shows the 
center of rotation 0, projection S of center of rotation O on scan PQ and 
projection S 1 of center of rotation 0 on scan P'Q 1 . FIGS. 15 and 16 also show 
the intersection T of these two scans. For two consecutive scans, distances ST 
and S'T, shown in FIG. 16, will be negligibly small (they are equal as follows 
from geometry) . Also, distance PS is equal to distance P'S' from the 
definition of 0 as the rotational center. Therefore distance PT can be assumed 
to be approximately equal to distance P'T. Since T is a point shared by two 
rotational transducer scans, a scan line passes through T perpendicular to the 
plane of the illustration and is shared by those two rotational transducer 
scans. Therefore, a comparison of the data in rwo sran lines S in two 
consecutive scans (vertical sran lines S in both srang shown in FIGS. 17 and 
18), will show one value S, 1 . ltoreq. S . ltoreq.M, for which these two data 
arrays will be most correlated (the presence of noise prevents these data 
arrays from being equal, even if they represent sran data for one and the same 
line in physical space) . Here any correlation function can be used. For two 
data arrays X=(x.sub.l, . . . , x.sub.n), Y=(y.sub.l, . . . , y.sub.n), the 
formula ##EQU1## 
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The graph of FIG. 19 shows the generic view of this correlation coefficient as 
a function of scan line number S. In this instance, X is the vertical scan 
line number S in the scan k and Y is the vertical scan line number S in the 
scan (k+1) , k=l to 178. For each number k, S(k,k+1) is computed as a scan line 
number which maximizes the oorrel ati on coefficient. For multiple maxima, there 
are several choices: a) take the minimum of those S maximizing correlation ; b) 
take the maximum of S; and c) take average of S. For example, if the total 
number of vertical scan lines in each scan were 11, then using method a) for a 
list of nnrrfilahinn coefficients {0.2 0.3 0.2 0.2 0.4 0.5 0.8 0.8 0.7 0.6 0.4} 
for scan lines 1 to 11 in scans 5 and 6, S(5,6)=7 is first computed. Then the 
average of all such S(k,k+1) is computed: ##EQU2## 

DEPR: 

To compute offset_y, a preferred embodiment utilizes the coherency between 
far-spaced scans. The doubly cross-hatched area in FIG. 14 represents area 
swept twice: once by the leftmost part PO" and once by the rightmost part 0"Q 
of the transducer array. FIG. 2 0 shows that the sought offset is equal to the 
distance O'O. For all the points T of the doubly cross-hatched area of FIG. 14 
lying at the same horizontal level as 0 ' , two slices PQ and P'Q 1 are found 
intersecting at T (see FIG. 21) . If slice PQ has scan number k, then slice 
P»Q' will have scan number (180-k) , since . angle . QTO ' = . angle . 0 1 TP ' . Knowing 
that P t O"=PO=&Scirc ;, then P'T=&:Scirc ;+0"T and PT=&Scirc ;-OT. The fact that 
0T=T0" leads to the following algorithm for computing offset_y. Defining 
slices number k and number (180-k) as conjugate slices, then for all conjugate 
pairs of slices, a correlation coefficient can be computed for the vertical 
scan line &Scirc ; -OT in slice k shown schematically in FIG. 22, and the 
vertical scan line &Scirc ; +OT in slice (180-k) shown schematically in FIG. 
23. A graph for this correlation coefficient as a function of OT for all 
possible values of OT can be plotted as shown in FIG. 24. Since the diagram of 
FIG. 14 represents one of four different combinations of signs of offset_x and 
offset_y, the value of offset OT (T in FIGS. 21 and 22) should be allowed to 
be both positive and negative, with the condition that it has an opposite sign 
in the conjugate slice. As previously with offset_x, the offset is selected 
that maximizes the correlation coefficient. Now, if offset OT corresponds to 
the maximum correlation coefficient, then for the pair of conjugate slices k 
and (180-k) 

CLPR: 

16. The method as recited in claim 15, wherein the step of employing the 
acquired scan data to calculate offset data comprises employing the acquired 
scan data to calculate a first offset along a first coordinate direction and a 
second offset along a second coordinate direction perpendicular to said first 
coordinate direction . 

CLPR: 

18. The method as recited in claim 17, wherein the step of employing the 
acquired scan data to calculate a second offset comprises the steps of: 

CLPV: 

means for calculating correlation coefficients for respective scan lines of 
each pair of rotational scans of said transducer array having scan numbers k 
and (k+1) , where k varies from 1 to K and K is the total number of scans minus 
1; 

CLPV: 

means for determining a scan line number S(k,k+1) corresponding to a maximum 
cor re 1 ati on coefficient for each scan number k; and 

CLPV: 

means for calculating correlation coefficients for a vertical scan line &Scirc 
; -OT in scan k and a vertical scan line &Scirc ; +0T in a scan conjugate to 
scan k for each pair of conjugate rotational transducer array scans, where 
line segment OT is defined by a point T where said pair of conjugate scans 
intersect and point 0 where a line 00 1 perpendicular to scan k and 
intersecting said center of rotation O' intersects scan k; and 

CLPV: 

means for determining a value for OT corresponding to a maximum cor re 1 ati on 
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coefficient for each said pair of conjugate scans; 
CLPV: 

wherein said second offset is calculated as a function of said value for OT 
corresponding to a maximum rnrrelahinn coefficient. 

CLPV: 

calculating correlation coefficients for respective scan lines of each pair of 
rotational transducer scans having scan numbers k and (k+1) , where k varies 
from 1 to K and K is the total number of scans minus 1; 

CLPV: 

determining a scan line number S(k,k+1) corresponding to a maximum rorrel at-i on 
coefficient for each scan number k; 

CLPV: 

calculating corr-eTati nn coefficients for a vertical scan line &Scirc ; -OT in 
scan k and a vertical scan line &Scirc ; +0T in a scan conjugate to scan k for 
each pair of conjugate rotational transducer array scans, where line segment 
OT is defined by a point T where said pair of conjugate scans intersect and 
point O where a line 00 1 perpendicular to scan k and intersecting said center 
of rotation O 1 intersects scan k; 

CLPV: 

determining a value for OT corresponding to a maximu m cor relation coefficient 
for each said pair of conjugate scans; and 

CLPV: 

calculating a second offset as a function of said value for OT corresponding 
to a maximum correlation coefficient. 



calculating cor re Tat I on coefficients for respective scan lines of each pair of 
rotational scans of said transducer array having scan numbers k and (k+1) , 
where k varies from 1 to K and K is the total number of scans minus 1; 

CLPV: 

determining a scan line number S(k,k+1) corresponding to a maximum correlation 
coefficient for each scan number k; 

CLPV: 

calculating correlation coefficients for a vertical scan line &Scirc ; -OT in 
scan k and vertical scan line fcScirc ;+0T in the scan conjugate to scan k for 
each pair of conjugate rotational transducer array scans, where line segment 
OT is defined by a point T where said pair of conjugate scans intersect and 
point 0 where a line OO' perpendicular to scan k and intersecting said center 
of rotation O' intersects scan k; 

CLPV: 

determining a value for OT corresponding to a max-imn m mrrplation coefficient 
for each said pair of conjugate scans; and 

CLPV: 

calculating a second offset as a function of said value for OT corresponding 
to a maximum correlation coefficient. 

CLPV: 

calculating correlation coefficients for respective scan lines of each pair of 
rotational scans of said transducer array having scan numbers k and (k+1) , 
where k varies from 1 to K and K is the total number of scans minus 1; 

CLPV: 

determining a scan line number S(k,k+1) corresponding to a maximum correlation 
coefficient for each scan number k; 



calculating correlation coefficients for a vertical scan line &Scirc ; -OT in 
scan k and a vertical scan line &Scirc ;+0T in a scan conjugate to scan k for 
each pair of conjugate rotational transducer array scans, where line segment 



CLPV: 



CLPV: 
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OT is defined by a point T where said pair of conjuga^P scans intersect and 
point 0 where a line 00 1 perpendicular to scan k and intersecting said center 
of rotation 0' intersects scan k; 



CLPV: 

determining a value for OT corresponding to a maximu m corral at -i on coefficient 
for each said pair of conjugate scans; and 

CLPV: 

calculating a second offset as a function of said value for OT corresponding 
to a maximum correlation coefficient. 
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TITLE: Time domain filtering for NMR phased array imaging 



ABPL: 

A method and apparatus for combining NMR response data of a sample from a 
plurality of closely spaced RF receiver coils of an nmr phased array in the 
time domain to form a composite NMR image wherein each of the RF receiver 
coils receives a different respective one of a plurality of NMR response 
signals, each of which is evoked from a portion of the sample within a field 
of a respective one of the receiver coils. The response signals are 
conditioned to develop a plurality of data point signals corresponding to the 
magnitude of each of the respective response signals from each of the receiver 
coils at successive time intervals. The data point signals are convolved by a 
time domain representation of a field map of the respective one of the 
receiver coils generating the corresponding one of the response signals. The 
convolved signals are rnmhinp.d on a time domain point-by-point basis to 
produce a time domain representation of the composite NMR image of the sample. 



BSPR: 

The present invention relates to nuclear magnptir resonance (NMR) imaging and, 

more particularly, to methods and apparatus for combining the simultaneously 
received data from a plurality of radio- frequency (RF) coils of an NMR phased 
array in the time, rather than image, domain to produce a composite image 
having high signal-to-noise ratio (SNR) throughout the image . 

BSPR: 

The term " NMR phased array" refers to apparatus, such as shown in Roemer et 
al. U.S. Pat. No. 4,871,969 (the disclosure of which is incorporated herein by 
reference) , wherein a plurality of closely-spaced RF coils is employed for 
simultaneously receiving different NMR response signals from associated 
portions of a sample (such as a patient in medical imaging) and combining the 
separate data from each coil to produce a single composite NMR image of the 
sample. By overlapping adjacent coils and connecting each coil to the input of 
an associated low- input -impedance preamplifier channel, the high SNR of a 
single surface coil can be maintained over f ields-of -view (FOV) characteristic 
of remote coils. 

BSPR: 

Currently, composite imagfis for nmr phased arrays are reconstructed in the 
image domain by combining the individual -image contributions on a weighted, 
point-by-point basis after first acquiring the complete NMR images for each 
separate coil. The reason for acquiring the separate images first is that the 
optimum set of weights needed to maximize SNR when combining the separate 
signals to produce the composite -image is a function of position, and so 
varies from point to point . While the phase shifters and transformers of the 
setup shown in FIG. 6 of Roemer et al . U.S. Pat. No. 4,825,162 can be adjusted 
to provide a composite imagg in the time domain having a high SNR at any 
particular point, different weighting must be applied for each point in order 
to obtain good sensitivity over the whole image . Thus, the conventional 
approach is to first separately acquire the different NMR image from each coil 
before combining the different individual i mages , on a point -by-point basis, 
to form the composite image . 

BSPR: 

nmr phased array imaging as described in the '162 patent, therefore, has the 
drawbacks of requiring large amounts of memory to store the separate coil 
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images before reconstruction and of necessitating long-time delays between 
acquisition of the last data point and onset of the first display of the 
reconstructed image . 



BSPR: 

It is desirable in nmr phased array imaging to be able to combine the data 
from the separate receiver coils as it is acquired on a time domain, rather 
than -image domain, basis without sacrificing SNR resol ution . Combining the 
data as acquired will reduce the total memory requirements of the system since 
only one combined data set would have to be stored and, because only the 
rnmhined data set will have to be transformed at the end of scanning, will 
also reduce the time between end-of-scan and first appearance of the composite 
image . 

BSPR: 

Among the several objects of the present invention will be noted the provision 
of a method and apparatus for forming a composite NMR image with high SNR 
throughout the image ; the provision of a method and apparatus for NMR 
spectroscopy and NMR imaging using data rnmhi nat-ion in real time; and the 
provision of a method and apparatus which overcomes the aforementioned and 
other disadvantages of the prior art. 

BSPR: 

In accordance with the invention, a method is provided for combining the 
simultaneously received different nmr response signals from a plurality of 
closely-spaced, overlapping RF receiver coils of an NMR phased array in the 
time domain, to form a composite image that has high SNR throughout the image . 
A filter scheme is utilized to develop a composite data set in the time 
domain, wherein each time point of the composite data is formed on the basis 
of contributions from previous data points and future data points. The data is 
passed through filter arrangements having one-, two- and three-dimensional 
filters before the signals are summed together. Each filter dimension 
corresponds to filtering in one of the time dimensions of k- space, i.e., the 
readout direction; the phase encode direction; and, in the case of 
three-dimensional imaging, the second phase encode direction. Filter 
coefficients are chosen to combine the data in a way that is simultaneously 
optimal for providing a high SNR at multiple points of the composite image . 
With more terms added to the filter, the SNR can be optimized over the entire 
image . 

DRPR: 

FIG. 1 (prior art) is a schematic view of an arrangement employed in the 
conventional method to combine the signals from the overlapping coils of an 
NMR phased array using image domain data processing techniques. 

DRPR: 

FIG. 4 is a filter function derived for the coil of FIG. 2 after Fourier 
transformation into -image space of the truncated representation of FIG. 3. 

DRPR: 

FIG. 5 is a schematic view of an arrangement employed in a method of combining 
the coil signals of an NMR phased array in the time domain using filters in 
accordance with the present invention. 

DRPR: 

FIG. 6 is a pictorial representation of truncated convolution in the 
arrangement of FIG. 5 for a two-dimensional single slice image . 

DRPR: 

FIGS. 10A, 10B, IOC and 10D are views of exemplary reconstructions in the 
image space (FIG. 10A) and time domain space (FIGS. 10B-10D) , respectively, 
showing the effect of filtering in accordance with the method of the 
invention. 

DEPR: 

FIGS. 1 and 5 show an NMR phased array 10, such as described in Roemer et al . 
U.S. Pat. No. 4,825,162, of a plurality of radio -frequency (RF) receiver coils 
12 (coils 1 through N.sub.c) defining an imaging volume for the NMR imaging of 
a sample, such as for the nmr medical diagnostic imaging of a human spine. The 
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separate surface coils 12 are identically configured and are arranged in 
closely- spaced relationship with overlapping f ields-of -view (FOV) , but with 
substantially no interaction between adjacent coils. The coils 12 are adapted 
as part of the NMR imaging process to simultaneously receive a different one 
of a plurality of nmr response signals each evoked from an associated portion 
of the sample enclosed in the imaging volume. As shown, each coil 12 has its 
own processing channel 14 including receiver circuitry 15 and an 
analog-to-digital converter 16. FIG. 1 is a schematic representation of the 
conventional data processing set-up for constructing a different NMR i mage for 
each channel 14 of a sample portion from the NMR response signals received by 
the associated coil 12 for that channel 14, and for subsequently combining the 
plurality of different images thus constructed, on a point -by-point basis, in 
the imagg domain, to produce a single final nmp imagp of all sample portions 
from which an nmr signal was received by any of the coils 12. FIG. 5 is a 
schematic representation of the corresponding set-up for performing the image 
reconstruction in the time domain utilizing the principles of the present 
invention . 



DEPR: 

As described in Roemer et al . U.S. Pat. No. 4,871,969, the optimal combination 
or weighting of signals from the individual coils 1-N.sub.c in the array 10 to 
achieve a high signal-to-noise ratio (SNR) is dependent on the location (x, y, 
z) of a particular volume element (voxel) . This is because the signal of each 
RF receiving coil C.sub.i is sensitive to nuclear spins in proportion to the 
field B.sub.i created by the coil, whereas the noise is "white noise" 
uniformly distributed over the image . Hence, the resultant SNR is a function 
of position. 

DEPR: 

Assume that I.sub.l (x,y) is the complex image obtained by reconstructing the 
data received from coil C.sub.i, and B.sub.i (x,y) is the RF magnetic field 
produced by coil C.sub.i. The real part of B is the x component (in magnet 
coordinates as opposed to the screen coordinates of the image) of the 
transverse RF magnetic field and the imaginary part of B is the y component of 
the field. If noise rnrrplaHnns are ignored (which will have little effect on 
i mage quality) and all coils 12 have approximately the same noise, the 
rnmh-i nation of separate i mages I.sub.i that optimizes the SNR in the composite 
imagg is given by ##EQU1## where I(x,y) is the composite image . 

DEPR: 

The complex imagg is really the product of the RF receiving coil magnetic 
field and the spin density S(x,y) given by 

DEPR: 

The magnitude of I(x,y) in equation (1) can be expressed as the sum of the 
products of the magnitudes of the i mage, and the magnetic field maps: ##EQU2## 

DEPR: 

Equation (3) gives a basic form usable in image reconstruction methods. 
rnmh-in-ing imagps using equation (3) is particularly convenient because the 
phase shifts of the individual receivers do not have to be known, and the 
image reconstruction programs do not have to carry the complex data. 

DEPR: 

FIG. 1 shows schematically the conventional process employed for combining the 
data in the image domain. The NMR signal from each coil 1-N.sub.c is sent 
through its associated channel 14 for processing by its own receiver 15, 
digitized by its own analog-to-digital (A/D) converter 16 and stored in 
digital form in its own assigned memory 18. After acquisition is complete, the 
data from each coil channel is separately subjected to transformation by 
processing means 19 and then mmhinpd point -by-point into a single composite 
image at summation means 2 0 in accordance with equation (3) . 

DEPR: 

To derive the time domain filtering method of the present invention (FIG. 5) , 
it was recognized that the mmhinpd image obtained using the image domain 
method is simply the Fourier transform of the original time dependent data. In 
accordance with equation (1), the optimal combination of images I(x,y) (i.e., 
that giving high SNR over the whole image) is obtained by multiplying each 
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separate coil image ^sub . i (x,y) by its correspondiri^^F coil magnetic field 
map profile B.sub.i (x,y) at 21, before summing the results at 20 (see FIG. 
1) . From linear system theory, however, it is known that convolution in the 
time domain is equivalent to multiplication in the spatial domain. Thus, for a 
single slice of multi-slice data the time domain representation of the 
composite imagp can be given by the two-dimensional convolution integral 
##EQU3## where N.sub.c is the number of coils, v.sub.i ( t . sub . r, t . sub . . PHI . ) 
is the time dependent nmr voltage signal measured on coil i, b.sub.i 
( . tau. . sub. r, . tau. . sub. . PHI . ) is the inverse Fourier transform of coil i's RF 
field map, A (t . sub . r , t . sub . . PHI . ) is the inverse Fourier transform of the 
composite data set and t.sub.r is the readout time for each phase encode time 
t .sub. .PHI. . 



DEPR: 

For a finite set of discrete samples, the inverse Fourier transform A(j,k) of 
the composite imagp is ##EQU4## where v.sub.i (j,k) is a matrix of NMR 
voltages measured on coil i and b.sub.i (j,k) is the discrete Fourier 
transform of the field map from coil i. The first and second arguments are the 
sample indices in the readout and phase encode directions, respectively, and 
N.sub.r and N. sub.. PHI. are the number of samples in the readout and phase 
encode directions, respectively. 

DEPR: 

At first glance, equation (5) which combines the data in the time domain does 
not appear to offer any computational advantages over equation (3) which 
combines the data in the image domain. According to equation (5) , the total 
number of operations required to obtain a single k- space data point of the 
composite i magp is proportional to the number of pixels N in the image , where 
N=N. sub. . PHI . . times .N. sub. r. Thus, the number of operations required to 
construct the entire k- space representation of the composite image scales is a 
factor of N.sup.2. So, for large values of N, rnmhin-ing the data in image 
space rather than in time space would appear to require far fewer 
computations. This is because the data from each coil channel is subjected to 
Fourier transformation at 19 (FIG. 1) before combining as a simple weighted 
sum at 20, and the number of operations for a Fast Fourier transform (FFT) 
scales as a factor of N log (N) rather than N.sup.2. 

DEPR: 

The number of computations required for the convolution can be greatly 
reduced, however, through the recognition that the RF field map 21 (shown in 
FIG. 1) are relatively slowly varying quantities across the image and can thus 
be suitably represented in abbreviated form for time domain processing 
purposes. It has been observed that the inverse Fourier transform of the field 
map is concentrated near the origin in the time domain (k- space) and thus the 
b.sub.i (l,m) terms in equation (5) can be truncated to a kernel containing 
relatively few terms. 

DEPR: 

By way of example, FIG. 2 shows the magnitude of a calculated sensitivity 
profile of a typical surface coil. The calculation is for a 40 cm FOV with a 
12 cm square loop RF receiving coil located in a plane perpendicular to the 
image. The main magnetic field is horizontal. The magnitude of its 
corresponding inverse Fourier transform is shown in FIG. 3, which is a contour 
plot of the k- space representative of the field map of FIG. 2. The constant 
contours are designated by arbitrary numerical units scaled to a maximum of 
1.0, with the maximum being at the origin. Only the center 31. times. 31 pixels 
of the magnitude are shown. Although the sensitivity profile (FIG. 2) occupies 
about 24% of the image FOV, the k-space representation (FIG. 3) has 
significant magnitude contribution only in the center few pixels. FIG. 4 shows 
the construction of a filter function profile corresponding to the sensitivity 
profile of FIG. 2, after truncating the filter coefficient of FIG. 3 by 
setting the magnitudes of the k space representation of FIG. 3 to zero outside 
the central 9. times. 9 pixel matrix of points, placing a Hamming window (see, 
R. W. Hamming, Digital Filters [Hall] pp. 102-105) around the data to avoid 
ringing, and then Fourier transforming the result into image space. A visual 
comparison of the derived filter function profile of FIG. 4 with the original 
profile of FIG. 2 shows little qualitative difference, except near the coil 
wires themselves. This indicates that a kernel of 9. times. 9 pixels is 
sufficient to give a good reconstruction. The error near the wires (located at 
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side lobes and the central re^^n 



the intersection of the side lobes and the central region of sensitivity) 
occurs because the RF magnetic field varies rapidly there and thus contains 
high spatial frequencies. Away from the coil, the RF field varies slowly and 
the 9. times. 9 filter kernel matches the profile well. 



DEPR: 

FIG. 5 shows a system and process in accordance with the invention for 
^nmhin-ing the separate coil data from an nmr array to obtain a composite image: 
with good SNR resolution in the time domain, using the above filtering 
technique. The front end of each coil channel 14' has its own receiver 15 and 
A/D circuitry 16 for receiving and digitizing the separately received signal, 
the same as for the corresponding channels 14 of the i mage, domain processing 
set-up of FIG. 1. But instead of storing the separate NMR images of each 
channel in a separate memory location 18 and Fourier transforming at 19 prior 
to summing, as done in the system of FIG. 1, the arrangement of FIG. 5, in 
accordance with the invention, filters the data with a field map filter 23 as 
it is acquired, and then sums the filtered data at summation means 20 1 prior 
to storing a pretransf ormation combined image in a memory 24. A single Fourier 
transformation is then undertaken by fast Fourier transformation means 25 to 
give the final composite image . The filters 23 provide the weighting necessary 
for summing the separate contributions from the channels 14 ' to give a good 
SNR rpsnhifinn in the end image. The filters 23 perform the operations defined 
by equation (5) . In contrast to the image domain data combination method 
employed by the system of FIG. 1, only one (or, possibly, two to obtain a 
uniform noise image.) Fourier transformation is required at the end of the 
scanning operation to produce the combined image. Thus the time domain 
filtering scheme of FIG. 5 avoids the large time delays from end-of-scan to 
first imagp appearance inherent in the process employed by the system of FIG. 
1. Moreover, the transformation process for the arrangement shown in FIG. 5 is 
independent of the number of coils utilized. Also, in contrast to the FIG. 1 
image* domain approach, the time domain method employed by the system of FIG. 5 
of the invention has the additional advantage that the data is combined in 
real time, as it is being collected, thereby reducing the data storage 
capacity necessary for each coil channel. 

DEPR: 

Various system architectures are possible for implementation of the time 
domain image reconstruction method of the invention. If the k- space 
representation (FIG. 3) of the field map for each coil 12 is reduced to a 
kernel size of N.sub.cr . times ,N. sub. c. PHI . in the readout and phase encode 
directions (see FIG. 4) , the k- space representation of the composite image may 
be given by ##EQU5## According to equation (6) , each data point from each coil 
12 contributes to a rectangular sub region of the composite k- space matrix. 
This relationship is shown pictorially in FIG. 6. 

DEPR: 

FIG. 6 is a pictorial representation of truncated convolution for a 
two-dimensional single slice image . The digitized signal 28 from a single coil 
12 has a number of data points 29, 30 of magnitude V.sub.n, V.sub.n+1, etc., 
corresponding to samplings of the analog voltage signal 28 taken at successive 
time intervals t.sub.n, t.sub.n+1, etc. The data from each coil 12 contributes 
to a rectangular region in the composite k- space (time domain) matrix. The 
size of the rectangular region is N. sub. c. PHI. . times . N. sub . cr . As each new 
data point enters, the subregion moves in the readout direction by one pixel. 
The signal 29 from time tn contributes to a rectangular subregion 31, and the 
signal 30 from time t.sub.n+1 contributes to a rectangular subregion 32 
shifted in the data memory 24 in the readout direction by one column. For each 
new phase encode step, the region moves down by one row. 

DEPR: 

A straightforward hardware embodiment of filter 23 for implementation of 
equation (6) is depicted in FIG. 7. The arrangement uses discrete components 
and needs a modest amount of memory to temporarily save the most recent 
N. sub. c. PHI. . times. N. sub. r data points from each slice, echo and coil. The 
multiplication element A performs the multiplication in equation (6) . The two 
innermost summations of equation (6) are done by summation element B and the 
temporary regi ste>r 36. The outermost sum is performed by summation element 
20». For each coil channel 14', a temporary storage memory 34 is connected to 
receive the digitized output signal of the associated coil 12. The memory 34 
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functions to save the data points 29, 3 0 for the successive time increments 
tn, tn+1, etc., until a sufficient amount of data is accumulated to complete 
one i,k point in the composite matrix. Assuming that k-space is covered in a 
linear fashion, this occurs during the N. sub. c. PHI. ' th phase encode step. As 
each subsequent data point enters from the left (viz., at the input end), the 
appropriate N.sub.cr . times ,N. sub. c. PHI . data is extracted from filter memory 
34 and is multiplied by the corresponding filter coefficients stored in a 
convolution kernel memory 35, summed by use of a temporary register 36, and 
then sent to an output port 37 as a weighted input to the composite summation 
means 20'. To avoid overflow of the input memory 34, a line of readout data 
must be completed and sent to the output port for each new line that enters. 
To accommodate mnl t-.-i -h~H r;e data, the memory must save the most recent 
N. sub. c. PHI. . times. N. sub. r data points from each slice and echo. 



DEPR: 

The architectures shown in FIGS . 7 and 8 do not exploit any special purpose 
filter chips that might be available. A number of companies make chips that 
compute one -dimensional convolutions very efficiently, and it may be desirable 
to utilize such chips to make a compact and efficient filter. FIG. 9 shows one 
way of breaking the two-dimensional filter 23 (shown in FIG. 5) into a series 
of one-dimensional filters. If this is done, the data rate into memory can be 
reduced by a factor of N . sub . c . PHI . , the number of convolution filter points 
in the phase encode direction. For each coil 12, a line of data in the readout 
direction is passed through a one -dimensional convolution filter 38 a total of 
N. sub. c. PHI. times. Each pass through the one-dimensional filter utilizes a 
different set of filter coefficients obtained from the filter memory and 
corresponds to the multiplication and innermost summation of equation (6) . The 
outermost sum of equation (6) , the sum over the coils, is produced by 
summation means 20*". Summation means 20'" also performs an addition for the 
read-modify- write operation which is equivalent to the second summation in 
equation (6) . The signals at the outputs 37" of the filters are summed by 
summation means 20'", and the results are added to the k-space storage memory 
24 for the composite images. For each readout line that enters from the left 
or input side, N. sub. c. PHI. readout lines in the composite stage memory 24 are 
modified. The same hardware can be used to combine three-dimensional data 
simply by changing the order in which data is read out of filter memory, and 
the order in which the results are added to the composite memory. 

DEPR: 

One proposed digital signal processing chip employs 32 multipliers that 
operate in parallel and can convolve a 16 bit-by- 16 point complex kernel with 
a 16 bit data path at a rate of 800 ns per complex I/Q pair. For multi -slice, 
data requiring a 16 point-by-16 point kernel and acquired at a 100% duty 
cycle, this corresponds to a maximum digitization rate of 12.8 .mu.sec. For 
conventional multi-slice imaging, 512 complex I/Q pairs can be typically 
acquired in 8 msec or 16 .mu.sec per point. Thus, it can be seen that a single 
filter chip may be sufficient for each coil channel 14 1 , with some spare for 
overhead such as reloading the filter coefficients. 

DEPR: 

For three-dimensional imaging, the incoming data rates are the same as for 
mnl t-.i -slice data, but another dimension of filtering is required. Since the 
data is acquired over many minutes with nearly 100% duty cycles, huge amounts 
of data are processed. It is therefore not practical to place temporary memory 
in front of each filter, and thus the data must be combined as fast as it 
enters. A three-dimensional image with 512 readout points and an 8 msec 
readout time would require 16 filter chips to keep up with the incoming data 
rate. A reduction in the filter kernel from 16 to 8 in the two phase encoding 
directions changes the data rate through the filter by a factor of four and 
this only four chips would be required. This may cause some degradation in the 
SNR (initial indications are not much) but the SNR will still be better than 
one could obtain without using the phased array. 

DEPR: 

A four-coil array was used to demonstrate the methods for single slice 
saggital imaging of the human spine. The array was made of 12 cm coils 
overlapping in a row in a manner similar to that shown in FIG. 4 of the '162 
patent. The four coils 12 were placed beneath the patient in a linear array 
running in the vertical direction. Each coil 12 had its own receiver and 
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digitizer. The image FQV was 40 cm with a composite ma^ix size of 

512 . times . 512 pixels. After the data was separately collected, it was rnmhined 

in the time domain with filter kernel sizes ranging from 1. times. 1 (a simple 

sum) to 9. times. 9. For comparison, the data was also combined -in th e -image 

domain. 

DEPR: 

Following the general procedure described in the '162 patent, the data from 
each coil was combined in the image domain by separately Fourier transforming 
the data from each coil 12, weighting each image by its corresponding field 
map image f and then summing the results. The image was then normalized into a 
uniform sensitivity image . FIG. 10A shows the resultant 512 . times . 512 saggital 
image of the spine. 

DEPR: 

To combine the data in the time domain utilizing the method of the present 
invention, the filter coefficients were determined by calculating the complex 
field map (magnitude and phase) for each coil over the full 512 . times . 512 
image matrix. Expressions for the field maps were obtained by integrating the 
Biot-Savart Law over a conductor placed at the centerline of the coil. The 
complex field map for each coil was then inverse Fourier transformed into 
k- space and truncated to the desired size. To avoid ringing in the image, the 
resultant filter coefficients were windowed in two dimensions with a Hamming 
window of the type described in R. W. Hamming, Digital Filters, supra. 

DEPR: 

To avoid interference patterns in the image , the raw data was corrected for 
the different phases and gains of the NMR receiving channels. To do this, a 
transmit loop approximately 2 cm in diameter was placed about 4 cm above each 
surface coil. The transmit loop was driven from the local oscillator used by 
the NMR receivers and this was phase-locked to the receiver. The resultant 
signal amplitude and phase measured at the output of the receivers were used 
to calibrate each channel. 



DEPR: 

Using equation (6), the truncated filter kernel was then convolved with the 
phase corrected NMR data. The results were then summed over the coils and 
Fourier transformed. To obtain a uniform sensitivity image, the filter kernel 
was convolved with the conjugate of the filter kernel and summed over the 
coils. This normalization map was also Fourier transformed and divided 
pixel -by-pixel into the image . 

DEPR: 

FIGS. 10B, 10C and 10D show the results of a simple sum, a 5. times. 5 point and 
9. times. 9 point kernel. The simple sum image (FIG. 10A) corresponds to an 
image obtained using a single large coil. As expected, the simple sum image 
had poor SNR (approximately two times lower than the point-to-point weighted 
sum image of FIG. 10A) and poor suppression of motion and wraparound 
artifacts. The 5. times. 5 kernel image showed significant improvement. The 
5. times. 5 image. (FIG. 10C) was developed by passing the time dependent data 
from each coil 12 through a two-dimensional filter with a 5. times. 5 kernel 
before summing. A faint wraparound artifact is visible at the top. The 
9. times. 9 image (FIG. 10D) was constructed similarly using a 9. times. 9 kernel. 
The 9. times. 9 result had almost the same quality as the composite image (FIG. 
10A) developed using the image domain techniques, except for minor differences 
near the edges of the image . 

DEPR: 

The differences of the edges between the images obtained using time domain 
(FIGS. 10B and 10D) and image domain (FIG. 10A) methods are due to the 
wraparound of the filter. Ideally, the filter corresponding to the coil at the 
bottom of the image should have no significant contribution at the top of the 
imag e . However, a 9. times. 9 filter can be made to roll off in only about l/9th 
of the image FOV and, thus, some wraparound is unavoidable. Combining the data 
in the image domain, however, allows one to filter to the nearest pixel or 
1/512 of the image FOV. To obtain exactly the same result in the time domain 
would required a 512 . times . 512 filter kernel. 



DEPR: 





In the above example, the filter coefficients were determined for the time 
domain filter hardware by calculating the RF magnetic field for each pixel in 
each slice for each coil. The results were then Fourier transformed, 
truncated, and then windowed. This method may not be fast enough, however, to 
be practical in the clinical environment. The filter coefficients are a 
function of the slice position and the operator of the NMR imaging system 
selects the slice orientation and the number of slices on a patent -by-patient 
basis. Within a minute or two after the selection, the nmr instrument needs to 
be ready to take data. Using the method of the example, to determine the 
coefficients for 30 slices of 512 . times . 512 pixel images using the four-coil 
array, 120 two-dimensional 512 . times . 512 complex inverse Fourier transforms 
would be required. If each transform takes 3 seconds on an array processor, 
this part of the computation would take 6 minutes. The calculation of the 
magnetic fields would probably add a few more minutes, thus creating a 
built-in delay of about 10 minutes, which might be unacceptable. A more rapid 
means of calculating the filter coefficients from the known positions or the 
RF receiving coils may be needed. 

DEPR: 

There are a couple of approaches to speeding up the process of filter 
calculation. One approach is to precalculate and store filter coefficients for 
common slice locations, but this might be too restrictive. Another approach 
involves a Fourier transform on a smaller grid. Since the resultant set of 
filter coefficients in k-space will be truncated, it is not necessary to 
sample the RF magnetic field at each and every pixel in the image before 
inverse Fourier transformation. The matrix size for coefficient calculations 
can thus be reduced from, say, 512 to perhaps 50 or fewer pixels. In the above 
30 slice calculation, this would decrease the time by a factor of 100. 

DEPR: 

A further approach involves the determination of the full , three -dimensional 
representation of a particular coil's RF magnetic field for a known position 
of the coil. This is done for each vector component, i.e., B.sub.x, B.sub.y, 
and B.sub.z, of the magnetic field. Each component of the RF magnetic field is 
then inverse Fourier transformed, truncated, and saved on disk for later use. 
Since rotations in real space are simple rotations in k-space and translations 
in real space are phase shifts in k-space, the set of filter coefficients can 
be derived for any coil location or orientation from this original stored set. 
For three-dimensional imaging, the field maps can be simply rotated and 
translated in k-space and then windowed to the desired size. For multi-slice 
data, the three-dimensional k-space data can be first rotated and then 
translated according to the slice position. The data can then be collapsed 
into two dimensions before windowing. Such methods involve relatively simple 
operations on small matrices, so they can be accomplished quickly. 

DEPR: 

In yet another approach, the NMMR image data itself can be used as the basis 
for determination of the filter coefficients. This is analogous to a 
sum-of -squares image approach and relies on the fact that the image itself is 
a measure of the coil sensitivity. By acquiring the center of k-space first, 
the raw data itself (actually its conjugate) becomes the filter coefficient 
for each channel . 



separate coil channels of an nmp phased array in the time domain using filters 
to produce a composite image having high SNR throughout the image . When 
compared with methods that romh-inp the data in the image domain, substantial 
reductions in the reconstruction time and the amount of memory required in the 
nmr imaging system are realized. In this way, systems using many coils can be 
made more practical . 

DEPL: 

where the * denotes the complex conjugate and C is an overall scale factor. 
The complex conjugate enters equation (2) because increasing angles of the RF 
magnetic field are defined to be positive in the direction of rotation of the 
nuclei. Greater angles of the RF magnetic field correspond to time delays 
(negative phase shifts) and thus the nmp signal is proportional to the complex 
conjugate of the RF magnetic field. 



DEPR: 

There has thus been described a method 
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DEPC: 

Combination in the Tmagp Domain 
DEPC: 

romhi nation in the Time Domain 
CLPR: 

I. A method for romhi ning nmr response data of a sample from a plurality of RF 
receiver coils of an NMR phased array in the time domain to form a composite 
NMR image, comprising the steps of: 

CLPR: 

8 . The method of claim 1 wherein the step of convolving includes the steps of 
obtaining an initial nmp imagp representation of the sample from each of the 
RF receiver coils and substituting the initial NMR image representation for 
the time domain representation of each receiver coil. 

CLPR: 

9. The method of claim 1 wherein the NMR composite image comprises a 
spectroscopic image of the sample. 

CLPR: 

10. A method for combining NMR response data from receiver coils of an NMR 
phased array in the time domain using filtering to form a composite NMR image 
having good overall SNR rpsnluhion, comprising the steps of: 

CLPR: 

II. Apparatus for receiving and nomhining NMR response data from a plurality 
of RF receiver coils of an NMR phased array to form a composite NMR image of a 
sample, comprising: 

CLPV: 

(a) receiving at each of the RF receiver coils a different one of a plurality 
of NMR response signals, each of the signals being evoked from a portion of 
the sample within a field of view of a respective one of the receiver coils; 

CLPV: 

(d) romhi ning the signals obtained by the step of convolving on a time domain 
point-by-point basis to produce a time domain representation of the composite 
NMR image of the sample. 

CLPV: 

(b) substantially simultaneously receiving at each one of the coils a 
different one of a plurality of nmr response signals, each evoked from a 
portion of the sample within the f ield-of -view of that coil; 

CLPV: 

(f) combining the multiplied temporarily stored signals of each coil, on a 
time domain point-by-point basis, to produce a composite NMR image of the 
sample . 

CLPV: 

a plurality of receiver circuits each connected to a respective one of the 
coils for receiving a corresponding NMR response signal, each response signal 
being evoked from a portion of the sample within the f ield-of -view of the 
response one of the coils; 

CLPV: 

means connected to said temporary storing means for combining said convolved 
series to provide a time domain representation of a composite NMR image of the 
sample . 
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DOCUMENT- IDENTIFIER: US 5568400 A 

TITLE: Multiplicative signal correction method and apparatus 



BSPR: 

In image analysis the fundamental variable is usually distance in one or two 
dimensions although the two-dimensional Fourier transform, also known as the 
Weiner transform, and the Weiner spectrum which express the information in a 
two-dimensional spatial frequency domain are also prevalent. Multivariate 
images, such as three color video signals and many satellite images where each 
picture element is characterized by a multichannel "spectrum" and also images 
constituting a time sequence of information, provide additional dimensionality 
in the data. Alternatively, in image analysis, the images can be summarized 
into histograms showing distributions of various picture elements, where the 
variable is then a vector of categories, each representing a class of picture 
elements, e.g. various gray levels of pixels, or contextual classes based on 
local image geometry. For multivariate images, the additional multichannel 
information may be included in the contextual classification. Time information 
can likewise be included in the definition of the categories in the variable. 
The above descriptions of two-way images also apply to three-way tomographic 
images, e.g. in MRI and X-ray tomography. 

BSPR: 

A prior approach to minimizing these errors has been to omit those portions of 
the spectrum having large variability from the data used in the regression. 
This approach is sometimes difficult to apply, because it may require many 
trials and operator judgments, and it is only partially successful at best. In 
a variation of this approach, the range of the spectral data included in the 
average spectrum used to determine the offset and slope coefficients is 
restricted to the vicinity of a strong isolated apftrhral feature, such as a 
solvent absorption band, thereby limiting the magnitude of the residuals and 
improving the accuracy of the correction. This variation has been applied to 
correction of the effects of scattering within the specimen in transmission 
spectroscopy. In many cases, however, there is no strong isolated band 
available for determination of the multiplicative correction. A related 
problem occurs in measuring one material through another with the pathlength 
through each material unknown and variable. In either case, better means are 
needed to accurately separate additive and multiplicative effects. 

BSPR: 

The present invention includes using any linear multivariate estimator to 
determine the correction coefficients such as multiple linear regression, 
generalized least squares, maximum likelihood regression, robust regression, 
estimated best linear predictor, partial least squares, principal component 
regression, Fourier regression, co-variance adjustment, or non-Euclidian 
distance measures. 



The present invention also comprises, as option B, generating modified 
reference spectra of the interfering components that contain only those 
portions of the original reference spectra that are orthogonal to, and 
therefore uncorrelated with, one or more reference analyte spectra. The 
coefficients generated for these orthogona 1 spectra are not influenced by the 
presence or magnitude of analyte information contained in the raw data even if 
the analyte spectrum is not included in the modeling or the coefficient 
estimator does not inherently orthogonalize the components, so they may be a 
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more correct representation of the magnitude of the spectral effects of the 
interfering components. These more accurate coefficients are then used to 
scale the original reference spectra prior to subtraction from the input data 
and the correction proceeds as in option A. This option reduces or eliminates 
the error in analyte spectral contribution that would otherwise be caused by 
subtracting an incorrect amount of a spectrum which contains some information 
equivalent to analyte information. 

DEPR: 

Like the conventional Martens, Jensen, and Geladi multiplicative scatter 
correction (MSC) , it seeks to normalize every input spectra to some reference 
"ideal" or "average" state by additive and multiplicative normalization, after 
having estimated the offset and slope parameters by some type of regression 
against some reference spectrum over some selected wavelength range, and this 
reference spectrum may be of the same kinds as those employed in MSC. However, 
it extends conventional MSC by explicitly modeling the effects of anticipated 
additive interferences and by optionally utilizing nonlinear modeling in 
deriving the additive and multiplicative normalization. This in turn improves 
the accuracy of the multiplicative correction, it allows removal of undesired 
interferants already at the multiplicative preprocessing stage, and it 
simplifies a causal understanding of the multiplicative correction and 
facilitates its interactive graphical optimization. It may also create 
interference reference spectra nrhhngnnal to the analyte (s) spectra for use in 
the modeling to avoid the effects of intercorrelation between the interferant 
spectra and the analyte spectra which would otherwise cause inaccuracies in 
estimating the interferant coefficients and in the subsequent subtraction of 
their contribution to the spectral data being normalized. 

DEPR: 

The R.sub.ka and R.sub.kj reference spectra are often statistical estimates 
extracted from the measured data from sets of specimens, although directly 
measured spectra are also useful in many cases. Honig ' s spectral 
reconstruction (D. E. Honigs, G. M. Hieftje, and T. Hirschfeld, A New Method 
for Obtaining Individual Component Spectra from Those of Complex Mixtures, 
Applied Spectroscopy, 38(3), pp. 317-322, a copy of which being annexed 
hereto) provides a method for extracting spectra from a set of mixture 
specimens based on knowledge of the concentration values. Principal component 
analysis (PCA) and partial least squares (PLS) provide orthogonal sets of 
spectra representative of the variation in the data. Stark's method (U.S. 
patent application Ser. No. 07/319,450) provides reference spectra for 
previously unknown variations based On analysis of replicate data. In the 
simplest operation of the present invention, i.e. correction for offset and 
multiplicative factors, the primary requirement of R.sub.ka and R.sub.kj is 
that they reasonably span the variation of X.sub.ki so as to stabilize the 
modeling and spectral accuracy and specificity are not essential. For the more 
complex options, in which R.sub.ka and R.sub.kj are incorporated into the 
output data as corrections, the quality of R.sub.ka and R.sub.kj become more 
important. The accuracy and specificity of the R.sub.ka spectral data is 
particularly important in orthogonalization of R.sub.kj either explicitly or 
implicitly and when used for added weight as described below. The spectral 
information in R.sub.ka and R.sub.kj may be represented in various ways with 
respect to redundancy and collinearity, for example one individual vector for 
each phenomenon, several replicates or specimens, statistical summaries 
(averages, bilinear components, square root of covariance matrices, etc.), or 
rotated representations where some or all collinearities have been eliminated. 
In a preferred embodiment, redundancy is eliminated by averaging so that the 
number of vectors equals the number of phenomena being modeled. 

DEPR: 

In a preferred addition to the basic preferred embodiment described above and 
shown in FIG. 3, additional correction spectra [R.sub.ka *b.sub.ai ] and 
[R.sub.ki *b.sub.ji ] are formed by matr-i x multiplier 340 from the reference 
spectra and their associated coefficients that were also generated in 
coefficient estimator 320. A matrix multiplier to form the spectrum [R.sub.kn 
*b.sub.ni ] for a single input spectrum i, consists of short term storage for 
the both inputs, a multiplication and summation circuit, an address sequencer 
which accesses the corresponding elements n of R.sub.kn and B.sub.nk and a 
second address sequencer which accesses the rows k of R.sub.kn and addresses 
the short term storage which keeps the resulting k. times . l_matxix. Mah-H y 
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multiplication is sFTso a standard function of available array processors. This 
additional combined correction spectrum may be used directly by subtractor 333 
to further correct Y.sub.ki, which becomes 



DEPR: 

The orthogonal component generator 360 provides for transformation of the 
reference spectra [R. sub. ks, R. sub.ka, R. sub. kj ] into a new set of spectra, 
[P. sub.ks, P. sub.ka, P. sub.kj ], some or all of which are> orthogonal to each 
other. If the reference spectra are latent variables derived from a single PCA 
or PLS analysis, they are orthogonal by definition. If they are measured 
spectra of components or otherwise separately derived, they will generally be 
intercorrelated, which if severe enough may cause errors in the coefficient 
values or failure of the coefficient estimator to complete its operation. If 
orthogonal reference spectra are created, new reference spectra may be added 
without requiring complete recalculation by the coefficient estimator. 
Orthogonal reference spectra also minimize the number of operations required 
_ by the coefficient estimator to determine the coefficients. 

DEPR: 

In a preferred embodiment, the nrt-hngnnal component generator performs a 
Gram- Schmidt orthogonal izat ion in accordance with ##EQU4## where I=the 
identity matrix 

DEPR: 

The variations of R. sub.ks are preserved in P. sub.ks, therefore the 
coefficient bs is not affected by the orthogonal izat ion. Each succeeding 
R.sub.kn is then orthogonalized against the matrix formed by the preceding 
orthogonal P.sub.kn spectra, until all reference spectra are orthogonalized 
into matrix P.sub.kn = [1, P. sub. ks, P. sub. ka, P. sub.kj ]. Each spectrum P.sub.kn 
comprises the residuals of the regression of R.sub.kn on the preceding 
orthogonal P.sub.kn. If a spectrum P.sub.kn is 0 or has only small values, it 
provides warning of dependence between spectra that could cause problems in 
coefficient estimation. In such case, the information is provided to the 
operator, or separate decision circuitry, to determine whether to delete the 
spectrum from the model, to downweight its importance, or to accept it without 
change. Orthogonal izat ion processing may be performed solely for the purpose 
of generating this warning information. When full orthogonal izat ion is chosen, 
the reference spectra input to mat-H x multiplier 340 are the P. sub.kj and 
P . sub . ka . 

DEPR: 

The orthogonal component generator and storage 360 comprises storage for 
P.sub.kn, the portion that is filled as the process proceeds comprising Z, 
storage for [Z'Z]-1, storage for the intermediate product Z(Z'Z)-1, storage 
for Z.sub.i, storage for the intermediate product Z'Z.sub.i, point by point 
multiply and sum logic, scalar inversion (1/a) logic, a subtractor, and the 
sequencer to select data from storage for processing, to control the 
processing sequence, and to direct storage of results. Circuit devices to 
perform these functions include the Intel 80287 math coprocessor for hardware 
implementation of the arithmetic functions, CMOS static ram chips (e.g. 4 
parallel Motorola MCM6226-30 128K. times . 8) to provide 32 bit resolution in 
storage of the digital data, and standard programmable array logic devices 
(PAL's) combined with a clock and counter as the sequencer. Each matrix 
element is acted on in sequence in accordance with the hardware logic. The 
required functions can also be obtained with a standard array processor 
operated in sequential fashion by the sequencer. 

DEPR: 

The contents of Z [Z ' Z] . sub . -1 is the transpose U' of U=[Z ? Z]-1Z' which is 
useful in finding multiple linear regression coefficients by matrix 
multiplication. 

DEPR: 

If even this degree of reference spectrum modification is undesireable , the 
orthogonal component generation is bypassed and the original reference spectra 
are passed to the coefficient estimator. 

DEPR: 

In the linear case, taking R.sub.kn - [1, R. sub . ks , R. sub . ka, R. sub . kj ], a matrix 
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i^j^ents observations at a value k^^^ 



where each row represents observations at a value k of the spectral variable 
and each column is a reference spectrum R.sub.kn incorporated in the model, 
coefficient estimator 320 fits X.sub.ki to Rnk by some method, minimizing the 
residuals e.sub.ki in 



DEPR: 

Methods for achieving this linear modeling include generalized least squares, 
maximum likelihood regression, robust regression, estimated best linear 
predictor, partial least squares, principal component regression, Fourier 
regression, mvarianrp adjustment, and others. For example, generalized least 
squares with generalized inverse models X.sub.ki by 

DEPR: 

A second preferred embodiment of the coefficient estimator which avoids matrix 
inversion is a principal components regression (PGR) device, which requires no 
pretreatment of R.sub.kn and no matrix inversion. 

DEPR: 

In the case of nonlinear modeling, the coefficient estimator 32 0 becomes more 
complex as each nonlinear coefficient becomes a vector of length k. C.sub.ki 
and D.sub.ki are therefore mahrir^R containing a number of coefficient vectors 
that depends on the form of the nonlinear model. 

DEPR: 

A preferred embodiment uses the coefficient estimator 320 illustrated in FIG. 
4 which employs Taylor series linearization. The model response generator 321 
calculates the vector F from the reference spectra Rkn, the present value Ari 
of the coefficients being generated by the iterative process, the variable k 
and the input spectral data Xki. This operation involves matrix multiplication 
and summation in accordance with the appropriate form of model as discussed 
above. The set of coefficients Ari, comprising cni of Cki, dni of Dki, and 
bni, are initially stored in coefficient Aq-1 storage 322a. They may be 
modified by means of adder 322b through addition of a weighted correction wGq 
or of increment dAr to one of the coefficients at a time to create the present 
values stored in coefficient Aq storage 322c and used by model response 
generator 321. The remaining functions will become obvious from the following 
description of the operation. 

DEPL: 

X is the spectral ordinate, e.g. absorbance, fluorescent energy, or pixel 
intensity or relative count, representing the measurement system response. The 
subscript i denotes the object or specimen while subscript k (k=l,2, . . . , K) 
is the spectral variable, k may be representative of a single dimension, e.g. 
wavelength in optical spectroscopy, two dimensions, e.g. time and wavelength 
in GC-IR measurements, or more depending on the measurement technology 
utilized. As used here, names of matrices are capitalized (e.g. R.sub.kj) and 
a matrix may comprise a single row Dr a single column of elements, for example 
X.sub.ki, Y.sub.ki, and R.sub.ks are individual spectra represented by single 
column matrices (vectors) of length K. Sets of multiple spectra form matrices 
(e.g. R.sub.ka, R.sub.kj, and Q. sub. km). Quantities that only exist as vectors 
or scalers are not capitalized. 

DEPL: 

Therefore, a general form for C.sub.ki and D.sub.ki can be described as a 
product of these series, ie a new series in terms of the powers of k and X and 
their cross products, removing redundant constants and terms in k or X and 
normalizing so that the linear magnitude information is kept in [R.sub.kn 
*t.sub.ni ] and C.sub.ki and D.sub.ki carry only the information relating to 
the nonlinearity . C.sub.ki and D.sub.ki are matri ppr containing k rows, and as 
many columns as required for the number of terms in the appropriate power 
series approximations. 

DEPL: 

Z comprises orthogonal columns therefore [Z*Z] is diagonal of size 

(i) .times. (i) and determining [Z'Z] -sup. -1 is trivial by inversion of the 

individual elements . 



DEPL: 

With this procedure, R.sub.ka may be omitted in the estimation of the b.sub.kj 



4 of 7 



12/03/2001 9:23 At 



coefficients without causing errors in their detemlnation. This method has 
the advantage of only removing analyte related information from the reference 
spectra, thus the analyte spectrum is unaffected, the interferant spectral 
shapes are minimally affected, and the coefficients have physical 
interpretations. In this case, the correct input to the matri x multiplier 340 
is R.sub.kj rather than P.sub.kj, to properly subtract the portion of R.sub.kj 
correlated with R.sub.ka. Implementation of this digital logic requires only a 
subset of the functions described previously. 

DEPL: 

where [ ] .about, means a generalized inverse and wh^re covari anrp m^friv v(i) 
can be iteratively updated based on the previous fit for this specimen i. 

DEPL: 

can be precomputed externally and stored with the reference spectra, thereby 
minimizing the requirements on the data normalizer 300. If full rank 
Gram-Schmidt orthogonalization is used, U.sub.nk is available from that 
process. In either case the calculation of the coefficients of the linear 
model involves a simple matrix multiplication. A matrix multiplier for 
U.sub.nk and X.sub.ki consists of short term storage for one or both inputs, a 
multiplication and summation circuit, an address sequencer which accesses the 
corresponding elements k of U.sub.nk and X.sub.ki and a second address 
sequencer which accesses the rows n of U.sub.nk and addresses the short term 
storage which keeps the resulting b.sub.ni values. In the more general case of 
multiple linear regression, matrix [R'R] must be formed and inverted to obtain 
[R'R] .sup.-l prior to matrix multiplication by R' to obtain U.sub.nk. This 
function can readily be accomplished with an available array processor and 
suitable logic sequencer. 

DEPL: 

where X=(x.sub.ik) is the matrix of spectral ordinates for i=l,2, . . . N 
objects, k=l,2, . . . ,K wavelengths, T=(t. sub.il) is the matrix of scores for 
objects i, bilinear factors 11=1,2, . . . ,L obtained from some bilinear model 
(e.g. principal component analysis, partial least squares, etc.), P=(p.sub.kl) 
are the loadings for objects i on bilinear factors 1, and E=(e.sub.ik) are the 
residuals between data X and model T*P' . The loadings P can then be decomposed 
into a function of a reference spectra r=(r.sub.k) (e.g. the Bean of the X 
data) and a matrix G=(g.sub.km) spanning the spectra for analyte and 
interference phenomena m=l,2, . . . ,M: 

DEPL: 

where d=(d.sub.l) and h=(h.sub.l) are vectors of length L, 1 is a vector of 
ones of length K, C=(c.sub.lm) is a matrix of regression coefficients of size 
L. times. M which quantifies the analyte and interference contributions, and 
F=(f.sub.lk) contains the residual loadings with the multiplicative, analyte, 
and interference phenomena removed, d, h, and C can be estimated by regression 
of P' on r, 1, and G by some method (e.g. weighted least squares) . C*G' could 
be reduced in size by elimination of effects if the relative size of the 
chemical or interferent effects are small 

DEPL: 

where U=is a matrix of eigenvalues and V* is a matrix of eigenvfirtors . By 
substitution, 

DEPL: 

The product T*U produces offset and interferent corrected scores and V is the 
matrix of corresponding spectra loadings associated with the corrected scores. 



DEPL: 

where S. sup.-l is the inverse of S. The fully corrected sr-nre matrix w is 
found in a similar fashion, 

DEPL: 

W=(W. sub.il) is the matrix of the offset, multiplicative, and interferent 
corrected scores which can be used as regressors in additive mixture models 
etc . 

DEPV: 
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Z=the matrix of vectors already transformed, 



DEPV: 

Z.sub.i T=transf ormed vector orthogonal to vectors already in Z. 
DEPV: 

I. Let X=(x.sub.ik) be the matrix of spectral ordinates for i=l,2, . . . N 
objects, k=l,2, . . . ,K wavelengths. The multiplicative effect is modeled 
from the spectral data using a standard multiplicative scatter correction 
(i.e., avoiding the use of components for practicing the present invention) 
yielding the corrected spectral data Z. 

DEPV: 

5 . Construct a new matrix of corrected spectra Z from X and repeat step 2 , 
step 3, and step 4 until convergence Occurs. 

DEPV: 

Let S be the diagonal matrix containing the elements of the product T*d. The 
fully corrected spectral data are found by 

CLPR: 

II. The method of claim 6 wherein said model includes a novariance adjustment 
technique . 

CLPR: 

18. The method as in claim 17 including the steps of generating modified 
reference spectra of interfering components that contain only those portions 
of original reference spectra of the interfering components that are 
nrthognnal to, and therefore uncorrelated with, at least one reference analyte 
spectrum, using the coefficients generated from said orthogonal reference 
spectrum to scale the original reference spectra prior to subtracting the 
scaled spectra from the data. 

CLPL: 

where X=(x.sub.ik) is the matrix of spectral ordinates for i=l,2 . . . N 
objects, k=l,2 . . . , K wavelengths, T=(t. sub.il) is the matrix of scores for 
objects i, 1=1,2 . . . L representing bilinear factors Obtained from a 
bilinear model, P=(P.sub.ki) are the loadings for objects i on bilinear 
factors 1, and E=(e.sub.ik) are the residuals between data X and model T*P» ; 

CLPL: 

where X=(x.sub.ik) is the matri x of spectral ordinates for i=l,2 . . . N 
objects, k=l,2 . . . , K wavelengths, T=(t. sub.il) is the matrix of scores for 
objects i, 1=1,2 . . . L representing bilinear factors obtained from a 
bilinear model, P=(P.sub.ki) are the loadings for objects i on bilinear 
factors 1, and E=(e.sub.ik) are the residuals between data X and model T*P ! ; 

CLPL: 

where X=(x.sub.ik) is the matrix of spectral ordinates for i=l,z . . . N 
objects, k=l,2 . . . , k wavelengths, T=(t. sub.il) is the matrix of scores for 
objects i, 1=1,2 . . . L representing bilinear factors obtained from a 
bilinear model, P=(P.sub.ki) are the loadings for objects i on bilinear 
factors 1 and E=(e.sub.ik) are the residuals between data X and model T*P'; 

CLPL: 

where X=(x.sub.ik) is the matrix of spectral ordinates for i=l,z . . . N 
objects, k=l,2 . . . , k wavelengths, T=(t. sub.il) is the matrix of scores for 
objects i, 1=1,2 . . . L representing bilinear factors obtained from a 
bilinear model, P=(P.sub.ki) are the loadings for objects i on bilinear 
factors 1 and E=(e.sub.ik) are the residuals between data X and model T*P'; 

CLPV: 

5) constructing a new matrix of corrected spectra Z from X; and 
CLPV: 

decomposing the loadings into a function of a reference spectral and, 
optionally, a matrix of spectral components for analyte and interference 
phenomena ; 
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CLPV: 

decomposing the loadings into a function of a reference spectrum and, 
optionally, a matrix of spectral components for analyte and interference 
phenomena ; 

CLPV: 

determining an S.sup.-l factor related to a matrix containing T*d for 
multiplicative effects from the new spectral scores T and using a factor U 
from calibration relating to additive and interferent effects to find 

CLPV: 

W=S.sup.-l *T*U; where W is the matrix of the additive, multiplicative and 
interferent corrected scores which can be used as regressors in additive 
mixture models; 

CLPV: 

fifth means for providing a signal representing the construction of a new 
matrix of corrected spectra Z from X; 

CLPV: 

means for providing a signal representing decomposed loadings which have been 
decomposed into a function of a reference spectra and, optionally, a matrix of 
spectral components for analyte and interference phenomena; 

CLPV: 

means for providing a signal representing decomposed loadings which have been 
decomposed into a function of a reference spectra and, optionally, a matrix of 
spectral components for analyte and interference phenomena; 
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